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Summary 

The postbuckling response of circular rings (or long cylinders) is examined. The rings are subjected to 
four types of external pressure loadings; each type of pressure is defined by its magnitude and direction at 
points on the buckled ring. Newton’s method is applied to the nonlinear differential equations of the exact 
inextensional theory for the ring problem. A zeroth approximation for the solution of the nonlinear equations, 
based on the mode shape corresponding to the first buckling pressure, is derived in closed form for each of 
the four types of pressure. The zeroth approximation is used to start the iteration cycle in Newton’s method 
to compute numerical solutions of the nonlinear equations. The zeroth approximations for the postbuckling 
pressure-deflection curves are compared with the converged solutions from Newton’s method and with similar 
results reported in the literature. 

Introduction 

The postbuckling behavior of circular rings (or of long cylinders) under external pressure is examined in 
this paper. The behavior is affected by the magnitude and direction of the pressure on the deformed ring. Four 
types of pressure loadings are considered: hydrostatic, constant directional, constant radial, and inverse-square 
radial. Each type is clearly defined in the paper in terms of the differential geometry of the deformed ring. 

The ring is modeled by exact inextensional shell theory. Newton’s method (ref. 1) is used to solve the 
nonlinear differential equations of the inextensional theory. Newton’s method leads directly to the classical 
linear buckling problem that determines the buckling pressures and their corresponding eigenfunctions. Next, 
by specifying the amplitude of the first eigenfunction, a zeroth approximation is provided to start the iteration 
in Newton’s method for the postbuckling solutions. A part of the zeroth approximation is to determine the 
pressure as a function of this amplitude by making the residual error in the zeroth approximation orthogonal 
to the first eigenfunction. This choice for the zeroth approximation in Newton’s method is essentially the same 
as the method of averaging, which is well-known in nonlinear vibration theory. Newton’s method is then used 
to correct the analytical zeroth approximation by a numerical solution of the linear variational equations of 
Newton’s method. The numerical solution uses numerical harmonic analysis to compute Fourier series solutions 
for the ring postbuckling deflections as a function of angular position on the undeformed ring. The accuracy 
of the numerical solution is nearly exact, since it is limited only by truncating the Fourier series and by small 
round-off errors. 

The theoretical results for the postbuckling of the ring are affected by terms in the theory that are often 
dropped in approximate nonlinear shell theory. Numerical results from inextensional theory are tabulated in 
the paper and compared with published results (refs. 2 through 10). 

The zeroth approximation from Newton’s method is compared with the numerical results for each type of 
pressure loading and proves to be a good approximation of the initial postbuckling behavior. The leading terms 
of the zeroth approximation are also compared with similar results in the literature derived from perturbation 
theory (refs. 2 through 6). The accuracy of the two approximations is comparable, but the analysis and 
algebraic operations for the zeroth approximation in Newton’s method are more direct than the analysis in the 
perturbation method for the ring problem. The direct calculation reduces the chance for algebraic errors that 
can go undetected in a perturbation analysis. 


Symbols 

A,C,D 

constants of integration in eigenfunction solution that depend 
on modal amplitude B 

a 

radius to point on undeformed circular ring 

B 

modal amplitude that appears as parameter in zeroth approx- 
imation, see equation (12) 

El 

bending stiffness of ring cross section 

is 

pure imaginary exponent appearing in eigenfunction 
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Bessel functions of first kind and order n 

stress resultants acting on cross-sectional area of ring 

nondimensional stress resultant, No* /El 

tangential and normal pressure, respectively, force/unit of arc 
length 

nondimensional components of pressure loading on deformed 
ring 

pressure magnitude, force /unit of arc length 
magnitude of position vector r 
position vector for point on deformed ring 
arc length of reference surface of ring 
unit vectors 

vector function that is solution of nonlinear differential 
equations of inextensional theory for circular ring 

vector function that is correction of U in Newton’s method 

magnitude of components of displacement vector, length 

nondimensional displacements, v = V/a and w = W/a 

maximum absolute value of w as function of 9 

coefficient of cos 29 in Fourier series solution for w\ parameter 
in results from perturbation theory 

angle between position, vector of point on deformed ring and 
corresponding vector for undeformed ring 

angle in cylindrical coordinates locating point on undeformed 
ring; independent variable in ring equations 

nondimensional load parameter, qa?/EI 

nondimensional variable, r/a 

angle between tangent vector at point on deformed ring and 
corresponding vector on undeformed ring 

coefficient of sin 29 in Fourier series solution for (f> 


Subscripts: 

H , R, C, S type of pressure loading: hydrostatic, constant directional, 

constant radial, and inverse-square radial, respectively 

j dummy subscript in Fourier series solutions 

m mth iteration cycle in Newton’s method 

Primes with symbols indicate derivatives with respect to independent variable 9. Superscript T indicates a 
transposed matrix. 


Geometric Relations for Types of Pressure Loadings 

The four types of pressure loadings considered in this paper for the ring, or for a long cylinder neglecting 
end effects, are defined in terms of the differential geometry of the deformed ring. The directions of the loadings 
acting at a point on the deformed ring are shown graphically in figure 1. 
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Figure 1. Pressure loadings applied to buckled ring. 


The position vector for a point on the neutral surface of the undeformed ring is the usual vector for polar 
coordinates, as indicated in figure 1: 


r = ae r 


(i) 


where 


e r = cos 9 i + sin 9 j 

and i and j are unit vectors along the x- and y- axes. The constant directional pressure loading is constant in 
magnitude and acts in the direction opposite to the radial unit vector e r , 


Q/2 ~ — <7# e r (2) 

The dimensions of the magnitude q R are force per unit of circumferential arc length* This pressure loading 
is normal to the undeformed ring, but does not necessarily remain normal at the corresponding point of 
application on the deformed ring. 

The other three types of pressure loading are initially applied normal to the undeformed ring, but the 
direction is a function of the deflection of the ring. The direction of the two radial pressure loadings is defined 
by the position vector for a point on the neutral surface of the deformed ring as follows: 

r = (a + W)e r + Ve 0 (3) 


where V and W are components of the displacement vector and e# is the unit tangent vector for the undeformed 
ring. The constant magnitude radial pressure loading acts in the opposite direction to the position vector r. 
The magnitude of the position vector is denoted by r, 

r = (r • r ) 1 / 2 = |(a + W) 2 + V 2 j ^ 

The unit vector for the position vector for the deformed shell makes an angle (3 with the position vector for 
the corresponding point on the undeformed shell, 


r 

- = cos 0 e r + sin fle# 
r 


(a + W)e r + Ve e 
r 


(4) 


3 



The constant radial pressure loading is defined by the relation 

q c = —q c ( cos/3 e r + sin/? e#) (5) 

The inverse square radial pressure loading is defined by the relation, 

2 

= -^j-(cos/? e r +sin/? e^) (6) 

and is the only loading considered whose magnitude varies as a function of ring displacement. 

The final pressure loading defined is hydrostatic pressure that acts normal to the deformed ring. The unit 
normal vector makes an angle </> with the undeformed position vector, 

n = — (cos *p e r -f* sin (p e^) 

The hydrostatic pressure loading is defined by the following equation: 


Q/Z — (7) 

The pressure loading terms appear in the nonlinear equilibrium equations for the ring or for the long cylinder. 
The notation for the stress resultants is shown in figure 2. The equilibrium equations are derived and 
nondimensionalized in appendix A. 



Figure 2. Stress resultants and pressures for buckled ring. 


Newton’s Method 

The nondimensional equations derived in appendix A (eqs. (A15)) are a sixth-order set of ordinary equations 
in the nondimensional variables n, </>, v, and w . The independent variable is 0, the angle in polar coordinates 
that locates points on the undeformed ring. The set of nonlinear equations is solved herein by Newton’s 
method following a general procedure outlined in reference 1. A solution of the equations where the ring 
remains circular can be written by inspection. This solution, the “membrane solution,” becomes the zeroth 
approximation in Newton’s method. A set of linear variational equations is derived in Newton’s method 
by expanding nonlinear terms in the nonlinear differential equations about the zeroth approximation and 
retaining linear terms. Periodic solutions of the linear variational equations exist for certain critical values 
of the nondimensional load factors for the pressure loading. The first critical value for each separate loading 
corresponds to the “buckling pressure” from classical linear buckling theory. 

After locating the critical load factors, Newton’s method is continued to seek postbuckled solutions near 
the buckling pressure. The continuation consists of fixing a zeroth approximation for the postbuckled solution. 
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The zeroth approximation has a prescribed amplitude that multiplies the mode shape determined by the 
periodic solution of the linear variational equations derived from the membrane solution. The new zeroth 
approximation is used to derive a new sequence of linear variational equations. The sequence converges to the 
desired postbuckled solutions. 

The convergence of the iteration for the postbuckled solution is improved by choosing a load factor in the 
new zeroth approximation that reduces the residual error in a least-squares sense. For the hydrostatic pressure 
case and the constant directional pressure case, the least-squares analysis can be carried out in closed form 
to give an approximate pressure-deflection relation. Similar accurate approximate results are obtained for the 
two radial pressure loadings, but the least-squares integrals cannot be evaluated in closed form. However, the 
integrals can be approximated by power series expansions or by numerical integration. The final results of 
the power series expansion are essentially the same as the results of the perturbation method reported in the 
literature (refs. 2 through 6); the difference in the two methods is the order of the mathematical operations. 

Linear Variational Equations 

The solution of the nonlinear differential equations is written as a vector function U in the following 
equation: 

u r = [n, <(>, V, u>] (8) 

Newton’s method solves the nonlinear problem by solving the recursive sequence of linear problems, with 
iteration counter m, as 

L(V m . 1 )6V m = -E(U m - 1) (9) 

Um = U m — 1 + 6\5 m (m = 1, 2, 3, • • •) (10) 

For the nonlinear ring problem (eqs. (A15)), the linear variational equations (eq. (9)), are as follows: 

f> n m + (! + 0m-l) ^0m + 0m- 1 ^0m + &Pt = ~{Et)m - 1 

(1 + 0m— 1 ) 6n m + n m-l <5 0m ~ + 6p n = ~(E n )m - i 

6w' m - 6v m + cos <£ m _i 6<j> m = ~(w' - v + sin 
6w m + 6v' m + sin 0 m -l <50m = ~{w + v’ - cos<£) m _i 

where 

(Et)m-i = W + (1 + 

(En)m — 1 ~ 4” <j> ) “ <t> 4~ Pn]m~ 1 

6p t = [Afl cos 0 + (Ac + \sp~ 2 ) cos (<p - /?)]m-l Hm 

+ [(Ac + Asp -2 )p -1 sin/3cos(0 - 0) - 2\ s p~ 3 cos/3sin((£ - 0)} m -\ 6w m 
+ [~(Ac + A 5 p _2 )p _1 cos/?cos(0 - 0) - 2Asp~ 3 sin/3sin(0 - 0)) m -\ 6v m 
s P n = -[A«sin0 + (A c + A s p' 2 )sin(0 - 0)) m -\ 6<f> m 

~ [(Ac + A SP~ 2 )P~ 1 sin/?sin(<^> - 0) + 2A s p~ 3 cos/3cos{<j) - 0)} m -\ 6w m 

+ [(Ac + A SP~ 2 )P~ 1 cos /? sin {<t> - 0) - 2\ s p~ 3 sin/?cos(<£ - 0)] m -\ Sv m 

All four nondimensional load parameters A#, A/j, Ac, and A $ appear in the linear variational equations. In 
the numerical analysis, any combination of the four pressures can be applied. The results in this paper are 
restricted to the four cases of each pressure acting alone. The linear variational equations are derived by 
expanding nonlinear terms in the Taylor series in the variables n, 4 >, v - and w and retaining linear terms in 
the corrections 6n, 8<j), 6v, and 6w. For example, 


(9a) 

(9b) 

(9c) 

(9d) 
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COS <f> = cos(0 TO _! + « cos^ m _i - sm4> m -i S<j> m 

P~ l ~ -Pm - 1 “ Prn—\ 6 P™ ~ ~Pm-l ~ Pm - 1 K 1 + w m-l) + «m-l <5^ m ] 
Another useful relation in deriving the linear variational equations is the differential of the angle /?, 

cos/3 6v — sin/? Sw 

dp = 

P 


Membrane Solution and Buckling Pressures 

Newton’s method requires a zeroth approximation to start the iteration sequence. The expansion about 
zero, 

= [ 0 , 0 , 0 , 0 ] 

gives the correction, 

(5Ui = [(5ni=-A//-A^-A(7“A5, 0, 0, 0] 


and the converged solution, U = £Ui, the membrane solution. The next step is to determine if other solutions 
of the nonlinear equations exist in addition to the membrane solution. The analysis in Newton’s method for 
identifying the existence of additional solutions consists of starting the usual iteration sequence with the known 
solution as the zeroth approximation. The linear variational equations in this case are homogeneous since the 
residual error of the known solution is zero. When the linear variational equations have nontrivial solutions 
that satisfy the boundary conditions, other solutions of the nonlinear problem exist. For the ring problem, the 
boundary conditions are that the solutions of the linear variational equations must be periodic in 9 for certain 
critical values of the load parameters. 

The solution of the linear homogeneous equations for the critical values of the pressure loadings is outlined 
in appendix B. The critical values of interest correspond to periodic solutions with two complete cycles around 
the ring, s = 2 in equations (B6): 


A# = 3 

(11a) 

Ar = 4 

(lib) 

>■0 = 1 

(11c) 

Xs = l 

(lid) 


Postbuckled Zeroth Approximation 

The critical pressure load factor in each case is a bifurcation point marking the intersection of an additional 
solution with the membrane solution. To compute these additional solutions using Newton’s method, it is 
necessary to start the iteration sequence with a new zeroth approximation that leads to a converged postbuckled 
solution. Selecting the new zeroth approximation for the ring problem proves to be a straightforward procedure. 
The analysis in appendix B determines eigenfunctions corresponding to each critical load. The eigenfunctions 
contain two arbitrary constants, and the zeroth approximation is determined by fixing these constants and 
adding the known function to the membrane solution. The new zeroth approximation results in a set of linear 
variational equations that no longer have periodic solutions for the homogeneous equations. On the other hand, 
the residual error on the right-hand side of the equations is periodic. Therefore, it is simple to show in theory 
that the particular solution of the equations is periodic and is a first correction for the zeroth approximation. 
An additional part of the analysis is to vary the load factor to keep the first correction small compared with 
the zeroth approximation. The choice of the load factor can be analytical or part of the numerical analysis 
of computing the first correction. The analytical choice will be considered first; the numerical analysis to 
determine the load factor is suggested by the analytical results. 

The load factor is determined by making the residual error vector orthogonal to the eigenfunction in the 
zeroth approximation. For the ring problem, this orthogonalization can be written explicitly since the geometric 
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nonlinear terms are primarily functions of the angle 0 . Therefore, guided by the form of the solutions with 
s — 2 in appendix B, let 

00 = B sin 20 (12) 

where B is now a prescribed constant rather than an arbitrary constant of integration and, as explained in 
appendix A and later in appendix B, an arbitrary phase angle in 0 is ignored. Once 0o 18 known, it can be 
substituted for 0 in the nonlinear equations. This substitution makes equations (A15c) and (A15d) linear 
equations in v and w. The particular solution of these equations for the zeroth approximation for v and w 
follows from the method of undetermined coefficients after making use of the expansions, 

00 

sin 0o = sin(B sin 29) = 2 ^ J2k-l{B) sin(2fc — 1) 29 (13a) 

*=l 


cos 0o = cos(i?sin20) = J${B) + 2^ J2j(B)cos4j9 

j = 1 

where Jk{B) represents Bessel functions of the first kind and integer order: 


(13b) 


w 0 = [-1 + MB)) + f; (2 f / 2fc-i (i?)cos(2fc - 1) 29 


00 2 

- Ify2 1 1 J 2 j( B ) cos 4 j9 (14a) 

oo 2 

w ° = S ( 2 k r !) ~ 2 _~x J 2k-l(B) sin(2fc - 1) 29 

oo g ■ 

- 2 E T^=- 7 J 2j(B)sinAj9 (14b) 

3-1 J 

The integration of the first nonlinear equilibrium equation (eq. (Al5a)) determines the zeroth approximation 
for the nondimensional stress resultant n as follows: 

= fp t d9 (15) 

where c is a constant of integration to be determined and the zero subscript has been dropped on 0. 

The only equation of the nonlinear equations not satisfied by the zeroth approximation is the equilibrium 
equation in the direction normal to the ring (eq. (A15b)). The residual error for the zeroth approximation for 
this equation is the result of substituting the expressions for 0 q and for and into the equation. The 

general expression for the residual is 


E n = (1+ 0') 


C — 0 / — 




A2 


- 0 W + 


Pn “ (1 + 0') j p t 


d9 


(16) 


The residual vanishes when 0 is part of an exact solution. For the zeroth approximation 0 q = B sin 20, the 
residual is an even periodic function of 0 and the fixed amplitude B. The constant of integration c and the 
pressure load factor A for any given loading type are free parameters. They are determined by equating to zero 
the first two terms in the Fourier cosine series for E n . Details of this orthogonalization procedure are contained 
in appendix C. The general form of the two equations that determine c and A is 
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(17a) 


c - 3J5 2 + Aao = 0 


£(2c + 6-3B 2 ) + Aai =0 


(17b) 


where ag and a\ are Fourier coefficients that depend on the type of pressure loading terms in the expressions 
for p n and pt and are functions of the modal amplitude B. Eliminating the constant c gives the load factor as 
a function of B : 


3(1 4- B 2 /2) 
a 0 “ d\/2B 


(18) 


Equation (18) takes the following forms for the four load types: 


1. Hydrostatic pressure: 



(19a) 


since ag = 1 and ai = 0 for this case. 

2. Constant directional pressure: 

A = 3(1 + g 2 /2) 

R Jq + BJ\(B) - J\(B)/2B 


(19b) 


which is derived from equation (13b) and the indefinite integral of equation (13a). For small values of B, 
the above equation is approximated by 




4(1 + B 2 /2) 
l + 3£ 2 /8 




(19c) 


3. Constant radial pressure: 

The Fourier analysis in appendix C does not give a simple expression for equation (18) for this case. The 
approximate results for the Fourier coefficients are 


_ . a\ 4 

“0 = 2and 2B = 3 


V B 2 

135 


so that equation (18) for constant radial pressure becomes 


x _ 3(1 + B 2 /2) 9 (, , 1W 2 \ 

C ~ (2/3) (1 + 16J? 2 /45) ^ 2 ( + 90 J 


(19d) 


4. Inverse-square radial pressure: 

The approximate results for the Fourier coefficients from appendix C for this case are 

= (I) s2 “ d S? = - (I) - (H) 

Equation (18) for inverse-square radial pressure is therefore 

_3(l_+fl»/2) /9\ / 9B2\ 

A5 ~(4/3)(l + 7 BV5) U)\ 10 J 


(19e) 
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This is the only loading with a negative coefficient for B 2 in the zeroth approximation for the postbuckling 
load; this implies that the ring is imperfection sensitive for an inverse- square radial pressure. 



Fixing the load factors completely defines the four zeroth approximations for solutions for the postbuckled 
ring. All the dependent variables of the problem are periodic functions of the independent variable 0 and 
contain the modal amplitude Basa parameter. 


Iteration for Postbuckled Solutions 

After determining the zeroth approximations, the residual error E n (eq. (16)) in the solution of the normal 
equilibrium equation is an even periodic function that has the coefficient of cos 40 as the leading coefficient 
in its Fourier series. The residual error can be reduced to zero as part of the usual iteration in Newton’s 
method that corrects the zeroth approximation. The iterative procedure cannot be carried out in closed form, 
but the insight gained in the analysis for the zeroth approximation makes clear the required properties of the 
numerical analysis. In fact, the algorithm is almost a direct numerical analog to the zeroth approximation. The 
linear variational equations for the ring (eqs. (9)) are solved for a periodic particular solution at each iteration 
step. For solutions at load factors near the critical values, the load factors are varied to keep these particular 
solutions for the corrections small. 

The nontrivial zeroth approximation derived in the preceding section is corrected by the iteration sequence 
in equations (9) and (10). The iteration sequence starts with the iteration counter m = 1 and 

Uq = [n o =no(B,0), </>o=i? sin 20, ^ 0 =^ 0 (£,0), wq=wq(B,6)} 

The load factor for each case is also a function of the constant B as derived in the preceding section. 

The linear variational equations (eqs. (9)) have variable coefficients and are solved numerically. The 
algorithm is a numerical analog to the zeroth approximation; the periodic zeroth approximation is corrected 
at each iteration step with a correction that is periodic with a least period of 7r in the independent variable 0. 
Formally, the iteration assumes that a periodic particular solution of the linear differential equations exists at 
each iteration step and that the corresponding homogeneous solutions can be ignored. 

Floquet theory applies to the solution of differential equations with periodic coefficients. This theory 
suggests that the algorithm used here will converge. Lack of convergence would be an indicator that the 
solution of the homogeneous equations cannot be ignored. However, the iteration converges rapidly for the 
practical range of load factors considered in this paper. 

The particular solution of the linear variational equations for the mth iteration step is sought in the form 


[SV ] m = 


r brim 
Swm 

L 6v m 


- J 

E n mj COS 2j0 
o 

J 

E <t>mj sin 2j0 
j = 1 

J 

E w mj COS 2jd 

j=0 

J 

E <i>mj sin 2.7*0 

- j=l 


(m = 1,2,3, •••) 


( 20 ) 


The form of the particular solution continues the form of the zeroth approximation where </> and v are 
odd functions of 0 and n and w are even functions. As discussed in appendixes A and B, the choice 
of even and odd periodic functions amounts to prescribing an arbitrary phase angle in the final periodic 
solution vector for the nonlinear system of equations. The series solutions are truncated to J terms 
with J as an input variable in the numerical solution. Substituting the assumed form of the particular 
solution (eqs. (20)) into the linear variational equations (eqs. (9)) and equating like coefficients of sin2j0 
and cos 2 j9 result in (4 J + 2) linear algebraic equations in as many unknown Fourier coefficients. In the 
numerical solution, the periodic coefficients in the set of differential equations are expanded in a Fourier 
series by using numerical harmonic analysis. These series in turn multiply the assumed series solutions 
for the dependent variables and the series for their derivatives. The products of series are combined 
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into one series for each linear variational equation. The Fourier coefficients of the combined series are 
linear combinations of the undetermined coefficients in 6 U (eqs. (20)), and finally these coefficients of like 
trigonometric terms are equated to known coefficients computed by expanding the residual error on the 
right-hand sides of the linear variational equations by numerical harmonic analysis. The final numerical 
operation is to solve the resulting set of linear algebraic equations for the unknown Fourier coefficients in 
£U. In summary, numerical harmonic analysis and trigonometric identities are used to implement a Galerkin 
solution. 

When the parameter B in the zeroth approximation is small, the algebraic equations to determine the 
Fourier coefficients in the particular solution are nearly singular. When B = 0, the equations are singular since 
they coincide with similar singular equations in appendix B. To keep the correction at each iteration step as 
small as possible, the load factor A is incremented in the residual error term for each type of pressure loading. 

The matrix form of the linear algebraic equations for the rath iteration step is then 


K m _iX m — E-m—i 6X 


'm - 1 


dE 

dX 


J m-1 


(ra = 1,2,3, •••) 


(21) 


where — \y m j , w m j , n m j, 

The algebraic equations (eq. (21)) are solved by a standard Gaussian elimination routine that factors the 
matrix K into an upper triangular form with pivoting about the maximum column element at each step of the 
factorization. The equations are factored so that the last diagonal element multiplies the correction on B, 
the coefficient of sin 29 in the series for <f>. This last diagonal element can be small, on the order of B, but the 
indicated division by the element is not performed. Instead, the increment on the load factor <5 A is computed 
so that the corresponding element on the right-hand side of the factored equations is zero. This operation is 
the numerical analog to the analysis in the zeroth approximation where B is fixed and A is determined as a 
function of B by an orthogonality relation (eq. (17b)). 

For larger values of B at loads away from the critical loads, the matrix K is well-conditioned, so that A can 
be fixed during the iteration and the value of B corrected during the iteration. 

Results 

Results for each type of pressure loading are presented in this section. For each case, the iteration to 
correct the zeroth approximation converged rapidly, including cases where the load factor A was prescribed on 
the postbuckling curve. The zeroth approximation for the pressure-deflection curves proved to be close to the 
final results. 

The results are also compared with other results in the literature. The zeroth approximations are compared 
with first approximations from perturbation theory. Also, the final results from the inextensional theory used 
herein are compared with results obtained by other theories. These latter results have significance in evaluating 
the accuracy of approximate shell theory for cylindrical shells under pressure loadings. 

The convergent pressure-deflection results for the four types of pressure loadings are summarized by the 
curves in figure 3. The nondimensional load factor A for each type of pressure loading is plotted as a function 
of the maximum radial deflection toward the center of the ring. Because of the sign convention, this deflection 
is denoted as |i^|max, the maximum absolute value of w as a function of 6 , and is determined by summing the 
Fourier series for w in the converged numerical solution. 

A summary of the accuracy of each zeroth approximation is plotted in figure 4. The dashed curves are the 
zeroth approximations for each type of pressure loading A as a function of the parameter B, the amplitude of 
the rotation (j) of the tangent vector in the postbuckled ring. The solid curves are the corresponding results 
for the pressure A as a function of B = 0i, the Fourier coefficient of sin 26 in the converged numerical solution 
for <f > . This comparison of results in the curves gives a quantitative measure of the validity of the assumption 
used to derive the zeroth approximation — namely, that the final postbuckled solutions are dominated by the 
eigenfunction corresponding to the buckling pressure. 

The accuracy of the zeroth approximation for the pressure-deflection relations suggests that design and 
optimization studies can be based on the zeroth approximation. Failure criteria can be based on stresses and 
strains in postbuckled rings or cylinders. The numerical solution shows that the zeroth approximation gives 
engineering accuracy for bending stresses and transverse shear stress resultants, but the terms in the Fourier 
series solutions with higher wave numbers affect the final numerical results for stresses more than the results 
for the deflections. 
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iype oi pressure 


Figure 3. Nondimensional pressure versus maximum radial deflection. 
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Figure 4. Zeroth approximation and converged solution from Newton’s method. 



Returning to the main results in this paper, postbuckled pressure-deflection relations, the results of the 
present analysis are compared with other results in the literature. Each case is considered separately in the 
following subsections. 

Hydrostatic Pressure 

Numerical results for this case are tabulated and compared with perturbation solutions given in table 1(a) 
and with an exact solution in table 1(b). The perturbation solution corresponds to the zeroth approximation, 
but a direct comparison of the two solutions is not possible. Formally, the solutions can be compared by equating 
the perturbation parameter w\ to 25/3, the coefficient of cos 29 in the eigensolution for w in appendix B. 
However, it is more consistent in comparing the load factor from perturbation theory with the final load factor 
from Newton’s method to set w\ equal to the coefficient of cos 29 in the converged Fourier series for w. This 
latter choice for w\ is used in computing all the tabulated results for the different pressure load factors from 
perturbation theory. 

The perturbation approximation for the hydrostatic case from Sills and Budiansky (ref. 2), El Naschie and 
El Nashai (ref. 3), and Budiansky (ref. 4) is 


81 

Atf = 3+32^? (22) 

The perturbation expansion tabulated in table 1(a) for the hydrostatic pressure case is slightly closer to the 
numerical solution for the load factor than the corresponding zeroth approximation from Newton’s method. 

For the hydrostatic pressure case, the normal and tangential equilibrium equations are not coupled with 
the displacements, and equation (15) is a nonlinear differential equation in the derivative of 0. Because the 
independent variable 9 does not appear explicitly in the equation, Carrier (ref. 5) was able to integrate the 
equation in closed form to give the load factor in terms of elliptic integrals. The results from Newton’s method 
are compared with Carrier’s results given in table 1(b). The load factors tabulated by Carrier are fixed input in 
Newton’s method and the maximum derivative of (j) as a function of 9 is compared from the two solutions. The 
small discrepancies in the results are due to truncation and roundoff errors in calculating the two solutions. 

An independent check of the numerical accuracy of the Galerkin solution used herein to solve the linear 
variational equations of Newton’s method is available for the hydrostatic case. The tangential stress resultant N 
must balance the pressure times the deformed radius at points on the ring where the transverse shear vanishes. 
This check on Newton’s method from statics shows at least an eight-digit accuracy for the results reported 
herein. 

The theory used herein and in reference 5 is classified as an exact inextensional theory. Long cylinders 
under external pressure can also be modeled with nonlinear shell theory. Rehfield (ref. 6) reports a decrease 
in hydrostatic pressure with increasing deflection for a perturbation solution of Sanders’ moderate bending 
theory (ref. 7). Hubka (ref. 8) checked Rehfield’s results by using Newton’s method. However, Hubka obtained 
pressure-deflection results closer to the inextensional theory by using a different approximate shell theory that 
accounts for the exact direction of the transverse shear stress resultant. 

Comparison of different theoretical models for predicting postbuckling response of long cylinders under 
hydrostatic pressure is beyond the scope of the present paper. However, the results reported here show that 
the question of which shell theory will yield practical results for cylinders under hydrostatic pressure is not 
closed. The theoretical models should also include the capability to determine “imperfection sensitivity” of 
cylinders with imperfections in shape. 

Constant Directional Pressure 

The numerical results for constant directional pressure are listed in table 2. The definition of the deflection 
parameters is the same as for the hydrostatic pressure case. The zeroth approximation in equation (19b) is 
very good for this case, much better than the result from perturbation theory of El Naschie and El Nashai 
(ref. 3). Rehfield (ref. 6) reports that shell theory predicts neutral stability with the pressure remaining constant 
with displacement in the postbuckling range. Singer and Babcock (ref. 9) point out that constant directional 
pressure is unstable dynamically for complete rings that are not prevented from rotating about the axis of 
symmetry. 
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Constant Radial Pressure 

Constant radial pressure results are summarized in table 3. The zeroth approximation from Newton’s 
method (eq. (19c)) and the perturbation result of Sills and Budiansky (ref. 2) 

X C = l + f 2 ™ l ( 23 ) 


show comparable accuracy. However, the perturbation results of El Naschie and El Nashai (ref. 3) are' not as 
accurate as the zeroth approximation. 

Inverse-Square Radial Pressure 

Results from radial pressure loading with the magnitude inversely proportional to the square of the radius 
are listed in table 4. This is the only loading that indicates imperfection sensitivity with postbuckled solutions 
at pressures below the critical bifurcation pressure. The zeroth approximation in equation (19d) is accurate, 
whereas the perturbation result of Sills and Budiansky (ref, 2) 


x 9 999 o 

As = 4 ‘ TH “■ 


(24) 


contains an algebraic error. The correct result (ref. 10) is 


. 9 999 2 

= 4 _ 224 Wl 


(25) 


One advantage of Newton’s method for the ring problem is that the analysis in the zeroth approximation can 
be checked for errors by using numerical harmonic analysis to determine the Fourier coefficients ao and a l 
that depend on the type of loading. In turn, the numerical harmonic analysis is part of the computation in 
solving the linear variational equations to correct the zeroth approximations. If there are programming errors 
in computing the residual error at each step of the iteration, the iteration to correct the zeroth approximation 
will not converge. Therefore, although Newton’s method is not foolproof, it is easier to detect errors in analysis 
or in computations than the perturbation method. 

Concluding Remarks 

The post buckling behavior of complete rings and of long cylinders under external pressure depends on the 
direction of the pressure loads as well as on the magnitude of the pressure. For rings and long cylinders, 
nonlinear inextensional theory predicts different buckling pressures and different postbuckling behavior for 
four types of pressure loadings that are all directed toward the center of the ring when the pressure is initially 
applied. Hydrostatic pressure that is applied normal to the deformed ring and constant radial pressure have 
stable postbuckling equilibrium configurations. Constant directional pressure is also stable if rigid body rotation 
of the complete ring is suppressed by an external support. The postbuckled equilibrium states for radial 
pressure with a magnitude proportional to the inverse square of the deformed radius are unstable and exhibit 
imperfection sensitivity. 

Newton’s method provides both analytical and numerical results for the inextensional theory. Analytical 
results for a zeroth approximation in Newton’s method compare favorably with perturbation results and can 
be derived more directly. Since much of the analysis can be written out explicitly, the ring problem is a good 
example of how the analysis in Newton’s method can guide in the proper choice of algorithms for numerical 
solutions that are both efficient and reliable. 

Results from inextensional theory presented herein and results reported in the literature from extensional 
shell theory with moderate rotations agree for the buckling pressures for long cylinders. The two theories show 
discrepancies in the postbuckling pressure-deflection curves for the case of hydrostatic pressure, the type of 
pressure usually found in engineering applications. 


NASA Langley Research Center 
Hampton, VA 23665-5225 
July 28, 1989 
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Appendix A 


Nonlinear Ring Equations 

The equilibrium equations for the ring are derived in this appendix and then written in nondimensional 
form. The equilibrium equations are written by summing forces in the two directions normal and tangential 
to the deformed ring. The tangential direction is determined by differentiating the position vector r for the 
deformed ring 

r = (a 4- W)e r + Ve$ (Al) 

with respect to arc length s. The unit tangent vector, the result of the differentiation, is 


t = 


~d,W dO' 

a 

ds ds 

e r i 


dV / rm d0 

^ + {a + w) Ts 




(A2) 


The tangent vector of the deformed ring makes an angle (p with the tangent vector of the undeformed ring 
e# so that, by definition, 

t = —sin >pe r + cos <pe# (A3) 

For any curve, the derivative of the unit tangent vector determines the curvature and its normal vector n, 


dt _ 
ds 7*1 


(A4) 


where 1/ri is the curvature. The derivative of the tangent vector for the deformed ring (eq. (A3)) is 


dt 

ds 


d<p dO 
ds ds 


^ (cos<pe r + sin (peg) 


(AS) 


The theory outlined here is inextensional. The differential of arc length ds for the deformed ring is assumed 
to be the same as for the undeformed ring, 

ds = a d9 (A6) 

Therefore the curvature for the deformed ring in the inextensional theory is 


d4 dh __ 1 (ty 1 

ds ds a dO a 


(A7) 


The additional relation 


dn _ 1 

ds 7*1 


completes the geometric relations necessary to derive the equilibrium equations for the ring. 


(AS) 


Equilibrium Equations 

The stress resultants for the ring are shown in figure 2. Summing forces on a differential element of the 
deformed ring leads to the vector equation 


d(N t) 
ds 


+ 


d(Qn ) 

ds 


4- P t t 4- P n n = 0 


(A9) 


where the tangential and normal components of the pressure loading are the projections of the pressure loadings 
defined earlier in the main text of the paper (also, see fig. 1): 


Pt = q R s\n<p 4- U c 4- j sin(<^ - /?) 


Pn=q H + Qr COS 4 + (*lc + ~3r') C0S (^ ~ P) 


(AlOa) 

(AlOb) 
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The scalar components of the vector equilibrium equation are 


Q + P* = 0 (Alla) 

as r\ 

-N + ^-+P n = 0 (Allb) 

r\ as 


The transverse shear force Q is eliminated by use of the moment equilibrium equation and the moment-curvature 
relation. The moment equilibrium equation is 



(A12) 


It 

in 


is assumed that the ring is thin so that the moment M is equal to the bending stiffness El times the change 
curvature of the ring. 



£7 M 
a d6 


(A13) 


Eliminating Q by using equation (A12) and then eliminating M by using equation (A13) give the equilibrium 
equations (All) as 


1 dN El 
a d6 a 3 



d 2 <t> 

W +Pt= ° 


1 

a 



AT El d 3 <t> „ „ 

N a 2 d03 + Pn ~ ° 


(Al4a) 

(Al4b) 


In addition to <f> and 0, the pressure loading terms l\ and P n are functions of the displacement components V 
and W . From the definition of the angle <t> in the unit tangent vector t in equation (A2) and the expression 
for t in equation (A3), 


sin ^ = “KS“'') < ai4c) 

“"■‘♦(t*;?) (A14d) 

The set of differential equations (eqs. (A14)) contains four equations in the four dependent variables 

and W in terms of the independent variable 6. Before solving the set of nonlinear differential equations, it is 

convenient to write them in nondimensional form. 


Nondimensional Equations 

Let the nondimensional variables be w = W/a, v = V/a, n = Na 2 /EI, p = r/a, and A = qa 2 /El. 
The nondimensional load parameters A will be given the same subscript as the dimensional factors q. The 
nondimensional nonlinear differential equations corresponding to equations (A14) are 


dn / d(b\ (fid 
le + \ 1 + de)de 2+Pt = 0 

(A15a) 

( 1 + ^)"“3^ +, ’» = 0 

(A15b) 

dw . , 

— — v + sin (f> = 0 

Clu 

(A15c) 
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(A15d) 


1 + u; + — - cos 0 = 0 


p t = X R sin 4> + + A SP 2 ) sin(4> - (5) 

(A15e) 

P n = A R COS <p + {\ C + A SP~ 2 ) COS {(f) -fi) + 

(Al5f) 

p 2 = (1 + w ) 2 + v 2 

(A15g) 

3 1 + w 

cos p = 

p 

(A15h) 

sin f3 = - 

(A15i) 


P 


The solutions of the nonlinear differential equations of physical interest are periodic in 6. If U = U(0) is a 
solution, then U(6 + 7) is also a solution where 7 is any constant phase angle. The value of the phase angle 
is arbitrary; knowing that the phase angle exists helps guide the application of Newton’s method by justifying 
the use of a sine series for certain functions in the numerical solutions and of a cosine series for other functions. 
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Appendix B 

Linear Variational Equations for Zeroth Approximation 

The linear homogeneous variational equations derived from the zeroth approximation, 


Uo = [no=-A//-A( 7 -A^-A 5 , 0, 0, 0] (Bl) 

are the following: 

6n j 4 Scj)^ 4 - (\r 4 A^ 4 \$) 6(f>\ — (A^ 4 A $) = 0 (B2a) 

bn\ — (Xff 4 4 A q 4 A 5 ) 6(f)\ — 6$" — 2 A 5 6w\ = 0 (B2b) 

6w[ — 6v 1 4 6cj) 1 = 0 (B2c) 

6wi 4 6v[ = 0 (B2d) 

The homogeneous equations have constant coefficients so that solutions can be sought in the exponential form, 

uf = [Ae a0 , Be* 6 , Ce aB , De af> j (B3) 

For a nontrivial solution, the constants satisfy the relations, 


C = -aD 

B = (a 2 + l) D 

A = jcr ^<r 2 + 1 ^ ^cr 2 + Xjj + Xr + Xq + A<j^ — 2 A 5 crj D 
and the exponents are roots of the characteristic equation, 

a 2 [a 2 4l) 4 o 2 {a 2 4 l) A# + (cr 2 + l) 4 cr 2 (cr 2 4 2 ) Xc 4 cr 4 As = 0 (B4) 

The characteristic equation is cubic in cr 2 and can be solved in terms of the pressure load factors. However, 
for periodic solutions of the homogeneous linear variational equations, the characteristic exponents must be of 
the form 


cr = ±is (z 2 = — 1 ; s = 0, 1, 2, 3, • ♦ • j (B5) 

Because of the condition that the solutions should be periodic, it is simpler to pick the characteristic 
exponents and solve the characteristic equation for the critical load factor. When each loading is treated 
separately, the critical load factors are 


A H = ( s2 ~ l) 

(B 6 a) 

A R = s 2 

(B 6 b) 

(s 2 - 1 ) 2 



(B 6 c) 

(* 2 - o 2 



(B 6 d) 
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The solutions of the differential equations for displacements corresponding to s = 0 and s = 1 are rigid 
body displacements and are not of physical interest, except for special cases where the ring has additional 
constraints. References 3 and 9 discuss this point for the case of A#, but the only solutions considered in the 
main text of this paper correspond to the exponents for s = 2. 

The bifurcation points defined by s = 2 are all isolated. This suggests that Newton’s method can be easily 
applied to computing solutions in the postbuckling range near the bifurcation points. The eigenfunctions 
corresponding to s = 2 have two arbitrary real constants when the exponential solutions with characteristic 
exponents a = ±i2 are put in real form. However, as noted at the end of appendix A, one of the two arbitrary 
constants is a phase angle in the periodic eigenfunctions. The other arbitrary constant is an amplitude of the 
periodic eigenfunctions and enters into the computation of solutions in the postbuckling range by Newton’s 
method. 
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Appendix C 


Fourier Analysis of Load Terms 

This appendix contains the details of the computation of the Fourier analysis of the load terms. The load 
terms fix the load parameter as a step in determining the zeroth approximation to the postbuckled solution 
for the ring equations. In addition to the relations between the angles <\> and /3 and the displacements v and w 
listed in the text of the paper, the two following relations are useful: 


sin (<f> - 0 ) = -p' 


COS (0 — /?) = p -f 


(1 4- w)v * — vw f 


Using the relations in equations (Cl), the load terms in equation (16) reduce to 

p n - (1 + <j>) J p t d9 = -f A r cos0 — (1 4- (f> ) J si n0 dd 

(1 4- w)v f — vw f 


4- A c 


4- A s 


2 p 4- p<t> f 4- 

P 

(1 + w)v f — vw l 

p 3 p . 


(Cla) 

(Clb) 


(C2) 


When </> is a periodic function of the angle 0 , the load terms are also periodic functions. The zeroth 
approximation for the load on the postbuckled ring is concerned with the two leading Fourier coefficients 
for each load term when 0 = 0 o = B sin 20. Squaring the zeroth approximation for (1 4-u>) and for v (eqs. (14)) 
and adding give 


P 2 - (1 4- w) 2 4- v 2 - 4- —Ji 4- 225 ^ ) 

+ ~ ^2 4 ^ cos 29 

4- — Jo«/i 4- - J\ 4- • ■ cos 40 4- • • ■ 

where the shorthand notation J n is used for the Bessel functions J n (B). Similarly, 

[(1 + w)v' - vw'} - (-fjf - + • ■ .) 

4- - JoJi 4- —J\J2 4- * cos 20 

4- cos 40 4- * * • 


(C3) 


(C4) 


The odd integer powers of p in the load terms are irrational functions of p 2 . These functions can be represented 
by power series expansions, but the resulting expressions are lengthy. In computing the load parameter 
corresponding to the zeroth approximation, retaining up through cubic terms in the constant B, the amplitude 
of 0q, provides a good approximation for continuing Newton’s method to compute postbuckling solutions. 
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The MACSYMA program (ref. 11) for symbolic manipulations was used to generate the following truncated 
expressions: 


r = i + — + 


rnZ m(m — 2 )Z 2 


+ m(m — 2)(m — 4 )Z 3 


(m = 1, -1, -3) 


where Z = p 2 — 1 and 


(C5) 


2 / 32 t2 \ / 32 . T 208 l3 32 _ _ \ 

Z - + ~ 27 J 1 “ 45 cos20 

+ (P J 1 3 ) cos4(? 

Z 3 = J l ) cos20 

where Ao = Jo ~ 1- The expansion for the coefficient of the load factor Ac is 

. (1 + w)v’ — vw' ( 4 14 o 8 o\ 

2 p + p<f> H = ^2Jq + ~ J 


(C6a) 

(C6b) 


( 2 


+ ( 2 BJ 0 + | J! - p 0 Ji + |a 0 Ji J 2 + 


+ T ~~ | JoJ i 3 - + S JlJ 2 ■ 1 J ° Jl Ja ) cos 20 (C7a) 

The leading terms in the power series that result from replacing the Bessel functions with their power series 
in their argument B are 


,, (1 + w)v f — vw f 

2 p + p4>* + = 2 + 


\3 135 J 


cos 2 9 


(C7b) 


The expansion for the coefficient of the load factor Ac; follows the same steps, 

- * - {>'■ - + + !**) 


(1 + w)v'-vw' 0' _ 8 2, 

“ Q “ T* I 


(-I-5M 


cos 29 


CO^h 

S?® 

+ 

CQ 

SI® 

1 


J 0 J 1 J 2 ) COS 20 

(C8a) 


(C8b) 
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Table 1. Summary of Different Results for Hydrostatic Pressure 


(a) Nondimensional pressure-deflection results 


Newton’s method 

Zeroth approximation 

Perturbation theory 



Rotation 


Deflection 



Deflection, 

coefficient. 


coefficient, 



| w | max 

B = 4>1 


w i 

Ah 6 

3.0 

0 

0 

3.0 

0 

3.0 

3.01128 

.06910 

.10 

3.01500 

.06659 

3.01122 

3.04543 

.14279 

.20 

3.06000 

.13275 

3.04461 

3.10346 

.22063 

.30 

3.13500 

.19804 

3.09928 

3.18710 

.30214 

.40 

3.24 

.26204 

3.17381 

3.30 

.38758 

.50085 

3.37628 

.32488 

3.26717 

3.40 

.45019 

.57275 

3.49206 

.36840 

3.34354 

3.44254 

.47424 

.60 

3.54000 

.38459 

3.37440 

3.62288 

.56381 

.70 

3.73500 

.44239 

3.49539 

3.84656 

.65504 

.80 

3.96000 

.49743 

3.62632 

4.00 

.70866 

.85816 

4.10466 

.52805 

3.70581 

4.12252 

.74743 

.90 

4.21500 

.54942 

3.76409 

4.20 

.77035 

.92468 

4.28255 

.56175 

3.79877 

4.40 

.82467 

.98304 

4.44955 

.59009 

3.88140 

4.46293 

.84043 

1.00000 

4.50000 

.59811 

3.90552 


“A# 


b *H 


3 1 + 



3 1 + 


27wl\ 
32 ) ' 


(b) Results for pressure versus maximum curvature change 



Rotation 
coefficient, 
B = <f>\ 

Maximum change of curvature, <//(0), for- 

Newton’s 

method 

Elliptic 

integral 

solution 

3.0 

0 

0 

0 

3.060 

.22949 

.47294 

.475 

3.227 

.43883 

.93239 

.959 

3.593 

.68487 

1.5175 

1.53 

4.076 

.88454 

2.0422 

2.13 

5.408 

1.19941 

3.0189 

2.98 
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Table 2. Summary of Different Results for Constant Directional Pressure 


Newton’s method 

Zeroth approximation 0 

Perturbation theory 

* R 

Deflection, 

1^1 max 

Rotation 
coefficient, 
B = <t> i 


Deflection 

coefficient, 

W 1 

A,r c 

4.0 

0 

0 

4.0 

0 

4.6 

4.0050 

.06926 

.1 

4.00501 

.06658 

4.01995 

4.02001 

.14337 

.2 

4.02011 

.13267 

4.07921 

4.04532 

.22184 

.3 

4.04556 

.19776 

4.17599 

4.10 

.34118 

.44377 

4.10112 

.28863 

4.37488 

4.12748 

.38964 

.5 

4.12927 

.32304 

4.46960 

4.18518 

.47781 

.6 

4.18881 

.38230 

4.65769 

4.25469 

.56804 

.7 

4.26126 

.43873 

4.86618 

4.33670 

.65971 

.8 

4.34766 

.49193 

5.08898 

4.40 

.72251 

.86796 

4.41496 

.52604 

5.24523 

4.50 

.81106 

.96349 

4.52229 

.57096 

5.46698 


“Tabulated values for A r for zeroth approximation are from equation (19b) rather than equation (19c). 
. 9w? 


Table 3. Summary of Different Results for Constant Radial Pressure 


Newton’s method 

Zeroth approximation 0 

Perturbation theory 

Ac 

Deflection, 

l^lmax 

Rotation 
coefficient, 
B = <p i 

Ac & 

Deflection 

coefficient, 

W\ 

Ac c 

4.5 

0 

0 

4.50 

0 

4.50 

4.50878 

.06936 

.10 

4.50657 

.06658 

4.50873 

4.53545 

.14377 

.20 

4.52717 

.13261 

4.53462 

4.58109 

.22269 

.30 

4.56447 

.19757 

4.57685 

4.64771 

.30555 

.40 

4.62312 

.26091 

4.63402 

4.73863 

.39176 

.50 

4.71021 

.32211 

4.70427 

5.00 

.56210 

.68949 

4.99615 

.43035 

4.86462 

5.23452 

.66427 

.80 

5.28623 

.48762 

4.96812 

5.50 

.74732 

.88899 

5.64187 

.52990 

5.05281 

6.00 

.84837 

.99701 

6.36480 

.57569 

5.15248 


a Tabulated values for for zeroth approximation are from power series in Bessel functions generated 
by the MACSYMA code and are slightly different from values from equation (19d). 


c\ 9 63 2 

Ac = 2 + w”?' 
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Table 4. Summary of Different Results for Inverse-Square Radial Pressure 


Newton’s method 

Zeroth approximation 0 

Perturbation theory 


Deflection, 
| w | max 

Rotation 
coefficient, 
B = 4 > i 

V 

Deflection 

coefficient, 

W\ 

V 

2.25 

0 

0 

2.25 

0 

2.25 

2.24 

.04870 

.07105 

2.23979 

.04734 

2.24 

2.17113 

.14350 

.20 

2.16969 

.13265 

2.17152 

2.07370 

.22213 

.30 

2.07117 

.19771 

2.07567 

1.93944 

.30464 

.40 

1.93656 

.26125 

1.94561 

1.77048 

.39048 

.50 

1.76884 

.32281 

1.78526 

1.56947 

.47909 

.60 

1.57136 

.38192 

1.59948 

1.33960 

.56989 

.70 

1.34770 

.43817 

1.39376 

1.08445 

.66233 

.80 

1.10133 

.49116 

1.17412 

.80785 

.75588 

.90 

.83519 

.54055 

.94687 

.50 

.85420 

1.0044 

.53799 

.58799 

.70810 


a Tabulated values for 
by the MACSYMA code 


\ s for zeroth approximation are from power series in Bessel functions generated 
and are slightly different from values from equation (19e). 
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